#梯度下降拟合ax+b
import numpy as np
import matplotlib.pyplot as plt

def real_func(x):
    return 2*x+1

def fit_func(p, x):
    f = np.poly1d(p)  # ps: numpy.poly1d([1,2,3])  生成  $1x^2+2x^1+3x^0$*
    return f(x)

def residuals_func(p, x, y):
    ret = fit_func(p, x) - y
    return ret


def dif(p,x,y):
    residuals=residuals_func(p,x,y)
    print(residuals)
    dif_alpha = -1*np.multiply(residuals,x).mean()
    dif_beta = -1*residuals.mean()
    return np.array([dif_alpha,dif_beta])

def grad(x,y,lamda=0.01,epoch=200):
    p_init = np.random.rand(2)
    X=np.insert(np.matrix(x).T,1,[1],1)
    for i in range (epoch):
        d_f = dif(p_init, x, y)
        print("epoch %d/%d d_f:"%(i,epoch),d_f)
        p_init=p_init+lamda*d_f
        output=fit_func(p_init,x)
        # while()
    return output


x_train=np.linspace(1,10,10)
y_train=real_func(x_train)+np.random.normal(0,0.5,x_train.shape[0])


# plt.title("")
# plt.plot(x_train,y_train,'g-*',label="real")
# plt.plot(x_train,grad(x_train,y_train,epoch=15),'r-*',label="SGD")
# plt.grid(True)
# plt.show()

'''多元线性回归'''
def real_func(x):
    return np.sin(2 * np.pi * x)

def fit_func(p, x):
    f = np.poly1d(p)  # ps: numpy.poly1d([1,2,3])  生成  $1x^2+2x^1+3x^0$*
    return f(x)

def residuals_func(p, x, y):
    ret = fit_func(p, x) - y
    return ret

def dif(p,X,y):
    #same
    # residuals=residuals_func(p,X[:,-2],y)
    residuals=X.dot(p)-y
    print("residuals",residuals)
    P=np.matrix(p)
    P=-residuals.dot(X)/len(residuals)
    print(P)
    return P

x_train=np.linspace(0,1,15)
y_train=real_func(x_train)+np.random.normal(0,0.1,x_train.shape[0])

def grad(x,y,lamda=0.4,epoch=150):
    M=5
    p_init = np.random.rand(M+1)
    X = x
    for i in range(2, M+1):
        X = np.column_stack((pow(x, i), X))
    X = np.insert(X, M, [1], 1)  # 多插入一列，用于+b
    print(X)
    for i in range (epoch):
        d_f = dif(p_init, X, y)
        print("epoch %d/%d\nd_f:"%(i,epoch),d_f)
        p_init=p_init+lamda*d_f
        print("p_init",p_init)
        output=fit_func(p_init,x)
        # while()
    return output

plt.title("")
plt.plot(x_train,y_train,'g-*',label="real")
plt.plot(x_train,grad(x_train,y_train,epoch=3000),'r-*',label="SGD")
plt.grid(True)
plt.show()

